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Abstract 

A concurrent multiscale modeling methodology that embeds a molecular 
dynamics (MD) region within a finite element (FEM) domain is used to study 
plastic processes at a crack tip in a single crystal of aluminum. The case of 
mode I loading is studied. A transition from deformation twinning to full 
dislocation emission from the crack tip is found when the crack plane is 
rotated around the [111] crystallographic axis. When the crack plane 
normal coincides with the [112] twinning direction, the crack propagates 
through a twinning mechanism. When the crack plane normal coincides with 
the [Oil] slip direction, the crack propagates through the emission of full 
dislocations. In intermediate orientations, a transition from full dislocation 
emission to twinning is found to occur with an increase in the stress intensity 
at the crack tip. This finding confirms the suggestion that the very high 
strain rates, inherently present in MD simulations, which produce higher 
stress intensities at the crack tip, over-predict the tendency for deformation 
twinning compared to experiments. The present study, therefore, aims to 
develop a more realistic and accurate predictive modeling of fracture 
processes. 
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I. Introduction 


Physics-based modeling of fracture begins at the nanometer-scale where atomistic simulation 
is used to predict the formation, propagation, and interaction of fundamental damage 
mechanisms. These mechanisms include dislocation formation and interaction, interstitial void 
formation, and atomic diffusion. The development of these damage mechanisms progress into 
microscale processes such as local plasticity and small crack formation. Ultimately, damage 
progression leads to macroscopic failure modes such as plastic yielding of components and large 
cracks exhibiting Mode I, II, and III deformation. 

Fracture inherently involves mechanisms that operate at a broad range of time and length 
scales - atomistic bond breaking at the crack tip, dislocation-governed plasticity in the damage 
zone near the crack tip, and statistically homogenized elasto-plastic properties of the material 
that influence the crack behavior at larger length scales. Multiscale modeling strategies have 
been developed to study these fundamental mechanisms while keeping the problem 
computationally tractable. Multiscale analyses attempt to bridge length scales by providing 
different physics-based models that can most appropriately represent damage mechanisms at 
each scale. These methods are especially useful in modeling of fracture in ductile metals, such as 
aluminum, which is the focus of this study. 

Recently, a number of atomistic simulation studies 1 ' 5 on crack propagation in aluminum have 
been published. The results of these investigations show that two main mechanisms of crack 
propagation operate. These mechanisms are propagation through deformation twinning, and 
propagation through the emission of full dislocations from the crack tip. However, some of these 
results disagree with experiment 2 . The experimentally observed dominant deformation mode is 
dislocation slip, while deformation twinning is rarely found. Whereas, the simulations 
significantly and systematically overestimate the degree of twinning in pure aluminum. This 
discrepancy between simulations and experiments has attracted considerable attention among 
researchers 3 ’ 5 because it prevents the reliable and accurate modeling of fracture in particular, and 
puts doubt on the reliability of the atomistic simulations in general. Most likely, the source of 
this discrepancy is related to the very different length (nanometers vs millimeters) and time 
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(nanoseconds vs seconds) scales at which simulations and experiments are usually performed. 
Nonetheless, the exact mechanism of how these length and time scales affect the propagation 
process remains unclear. 

In an attempt to resolve the issue, Tadmor and Hai 3 * systematically studied the role of 
crystallographic orientation of the crack plane on the propagation mechanism. Their assumption 
was that the constraint imposed by the boundary conditions (dimensionality constraints in a 
periodic system or surface effects in very thin specimens) 3 used for systems of very small sizes 
in atomistic simulations inhibit dislocation slip and favor twinning. 

More recently, Warner et al. 5 suggested that the crack propagation mechanism is strongly 
affected by the rate of propagation through the required applied stress intensity at the crack tip. 
To obtain a measurable propagation in the time frame of the simulation (a few nanoseconds), 
much larger stress intensities have to be applied as compared to the experimental time frame 
(seconds). According to Warner et al. 5 , higher stress intensities favor twinning, and thus 
twinning is more often found in simulations than in experiments. 

The present study investigates in detail the conditions responsible for twinning or dislocation 
emission from a crack tip under Mode I loading. A recently developed embedded statistical 
coupling method (ESCM) 6 ’ 7 for concurrent multiscale modeling is used. The approach is based 
on the construction of a coupled atomistic-continuum model. In this method, an atomistic 
domain is embedded within a larger continuum domain (Fig. 1). The circular atomistic region 
containing the crack tip is simulated using molecular dynamics (MD), while the surrounding 
continuum region is simulated using the finite element (FE) method. The authors have previously 
used the ESCM approach to study intergranular fracture along a grain boundary in aluminum 8 . 
The same approach is applied here to study transgranular fracture in a single crystal of 
aluminum. 
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Figure 1: Model geometry of the MD-FEM coupled system with an embedded edge crack ending inside the 
MD domain. 


II. The Simulation Approach 

The simulation approach used in this study, based on the ESCM, is a coupled atomistic- 
continuum model that allows a large atomistic domain (containing 10 4 5 or more atoms) to be 
embedded within a continuum domain of micron dimensions. The atomistic domain is simulated 
using MD, where the interatomic forces are represented by an atomistic potential suitably fitted 
to reproduce the material properties of aluminum 9 . The continuum domain is simulated by using 
the FEM with anisotropic elastic properties derived from the aluminum potential at room 
temperature (300K). The ESCM approach provides elastic boundary conditions for the atomistic 
domain by transferring the applied far field mechanical load through the continuum domain into 
the atomistic system. 
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The ESCM method is based on the concept that the continuum representation of a material, 
in terms of stress-strain fields, is a statistical representation of its atomic structure. Following this 
idea, the connection between the atomistic and the continuum representations at the MD/FE 
interface is performed using a statistical mechanics approach. The method uses statistical 
averaging over both time and volume of atomistic subdomains at the MD/FE interface to provide 
nodal displacement boundary conditions to the continuum FE domain. The FEM generates 
interface reaction forces that are uniformly distributed over the interface atoms in the form of 
constant traction boundary conditions 10,11 to the MD domain. Thus, this approach is based on 
solving the special boundary value problem (BVP) at the MD/FE interface for a coupled MD- 
FEM system and may be described as a local-nonlocal BVP because it relates local continuum 
nodal quantities with statistical averages of nonlocal atomistic quantities over selected atomic 
volumes. For this study, one finite element at the interface encompasses a region of several 
hundred atoms, the positions of which have been averaged over a 1 ps time interval. At this 
scale, the discreteness of the atomic structure is sufficiently homogenized so that the FEM 
domain responds to the atomistic domain as an extension of the continuum. The constant traction 
boundary conditions of the MD domain ensure that the elastic field from the FEM domain is 
correctly transferred to the atomistic region. An iterative procedure using the MD statistical 
displacements establishes a balance between the FEM-computed forces and the MD-computed 
forces at the interface. This force balance ensures stress continuity at the interface. 

The model geometry used to study transgranular crack propagation in an aluminum single 
crystal is shown in Fig. 1. A circular MD domain with a diameter */md = 40 nm was embedded in 
a square FEM mesh with side dimension dm = 20<4id = 0.8 pm. In this model, the system 
represents a single crystal with a pre-inserted edge crack propagating along the v-d beet ion. The 
crack plane normal is along the y-direction, and the crack front is along the ^-direction. The crack 
from the continuum region is initially extended into the MD domain equal to 1/3 c/md (see the 
enlarged central zone in the inset in Fig. 1). The atomically sharp crack tip in the MD domain is 
formed by screening the atomic forces between the atoms on both sides of the crack plane 
starting from the MD/FEM interface along a distance of 1/3 Jmd inside the MD domain. 
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The crack propagation direction, x, is fixed along the [ill] axis, so that the crack plane is 
perpendicular to the (l 1 l) slip plane. The crack front lies in the intersection of the (l 1 l) slip 

plane and the crack plane (Fig. 2). The orientation of the crack front line is chosen as the z- 
direction in the model. Under these circumstances, the only variable parameter for the 
crystallographic orientation of the crack is the twist angle, cp, formed between the crack plane 
normal (the y-dircction in Fig. 2) and the [112] twin axis, lying in the (ill) slip plane. The 

reason behind this crystallographic setup is based on the Peierls criterion for deformation 
twinning suggested by Tadmor and Hai 3 . According to this criterion, the tilt angle, 0, between 
the slip plane and the crack plane does not affect the tendency of the crack tip to nucleate a twin 
or a full dislocation under mode I loading. This tendency depends only on the twist angle, qp. In 
this way, setting 0 = 90° is a convenient choice, because it permits the cp angle to vary while 

keeping the crack front parallel to the most 
active (ill) slip plane. Because of the 

symmetry of the f.c.c. lattice, one can only 
examine variations of cp between 0° and 30°. 
Altogether six values of cp have been considered 
here: cp = 0°, 13.90°, 16.10°, 19.11°, 21.05°, and 
30° that satisfy the crystallographies of the 
system. The corresponding crystallographic 
orientations of the crystal lattice are presented in 
Table 1. These orientations were chosen in order 
to examine values for cp around 15°, as will be 
discussed in Section III. 


Figure 2: Crystallographic orientation of the crack with 
respect to the (ill) slip plane in the f.c.c. lattice. 

Periodic boundary conditions were applied in the z-direction of the MD system to emulate 
bulk atomic state. To achieve a correct calculation of the interatomic potential energy and force, 
the thickness h in the z-direction of the MD domain has to include a whole number of periods of 
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the crystal lattice in this direction. In addition, to ensure that at least 10 atomic planes are present 
in the z-direction (for consistency with the previous simulations by this group 4 ), a condition h > 2 
nm was imposed. The choices of h at 0 K for each value of cp are given in Table 1. The condition 
that the thickness in the z-direction is the same everywhere in the MD domain corresponds to a 
plane strain condition in the continuum domain. 
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Table 1: Crystallographic orientations of the crystal lattice and thickness, h, of the MD system for the applied set of 
cp - angles. 


The system has been loaded in uniaxial tension along the y-axis, through applying uniform 
displacement boundary conditions on the outer boundaries of the FEM square domain. The 
resulting stress intensity at the crack tip follows the continuum definition for K\ for the presented 
geometry 12 . The loading started from zero strain, 8 = 0, and was incremented by 0.1% every 20 
ps of simulation time, equivalent to a strain rate of 5x1 0 7 s' 1 . The applied strain increments of 
0.1% resulted in increments of the stress intensity at the initial crack tip of K\ = 0.049 MPa-rn . 
The strain rate was sustained until a nano-twin or a full dislocation was emitted from the crack 
tip. The average stress along the periodic z-direction in the MD domain was maintained at zero 
by applying the Parrinello-Rahman constant stress technique 13 . The MD system was simulated at 
constant temperature (300K) using the Nose-Hoover thermostat 14 . 

To analyze the processes taking place in the atomistic domain, the common-neighbor- 
analysis (CNA) technique was applied 15 . CNA identifies the local crystallographic state of an 
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atom based on the relative positions of its neighbors. For this study, two specific crystallographic 
states were of interest: the f.c.c. state - the ground state - and the h.c.p. state, which indicates a 
twin boundary (a single plane of h.c.p. atoms) or a stacking fault (a double plane of h.c.p. atoms). 
Atoms that were not identified as being in either of these two states were considered as 
disordered atoms. In a single crystal, the disordered atoms indicate a dislocation core or a 
vacancy. In addition, atoms that have lost more than 1/3 of its neighbors in the vicinity of their 
interaction range were considered as free surface atoms. In this model surface atoms indicate the 
crack surface. In the simulation snapshots presented in the paper (e.g., Figs. 3, 6 and 7, discussed 
in Section III), atoms in these various crystallographic states were color-coded to aid in 
visualization. All atoms in an f.c.c. state were colored in grey, and f.c.c. atoms that interface with 
the FEM nodes 6 ' 7 were colored in dark-gray. Atoms in h.c.p. state were colored in red, disordered 
atoms were colored in blue and surface atoms were colored in brown. For example, in Section 
III, the consecutive emission of the three partial dislocations, formations dl, d2, and d3 seen as 
blue dots in Fig. 3(a), produce a nano-twin indicated by parallel red lines that begin at the crack 
tip. 


III. Results and Analysis 

III. 1 . Twin-dislocation transition with orientation 

The first set of simulations was performed to study the initiation process of a twin or a full 
dislocation emission from the crack tip during crack propagation under mode I loading. As 
explained in Section II, the crack was gradually loaded by increasing the applied tensile strain in 
the FEM domain and monitoring the processes at the crack tip in the MD domain. The initial 
twin or dislocation nucleation for all six simulated crystallographic orientations is given in Fig. 
3. The nucleation process was first observed when the stress intensity reached K\ = 0.34 MPa- 
m . However, because of the incremental increase of K\ in steps of 0.049 (see Sec. II), the 
critical stress intensity K\ c, at which the first dislocation dl in Fig. 3 is nucleated, was estimated 
to be between 0.29 and 0.34 MPa-m 12 . 
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Figure 3: Simulation snapshots of the initiation of plastic processes at the crack tip for the set of 
orientations given in Table 1. 


To distinguish between the occurrences of twinning or full dislocation emission, at least two 
of the dislocations nucleated from the crack tip (dl and d2 in Fig. 3) need to be examined. If the 
second dislocation leaves an extrinsic stacking fault 15 behind, indicated by two parallel lines of 
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h.c.p. atoms separated by one line of f.c.c. atoms 16 (Fig. 3a and Fig. 3b), then a nano-twin 
nucleation has taken place 15 ' 16 . Further emissions of dislocations would most likely broaden the 
nano-twin 15 as the third dislocation has done in Fig. 3(a) and Fig. 3(b). If the second dislocation 
terminates the stacking fault left behind the first dislocation (Figs. 3c-f) then a full dislocation 
(consisted of two partial dislocations, d 1 + d2 in Fig. 3c) emission has taken place. 

The snapshots in Fig. 3 indicate a transition from twin nucleation to full dislocation 
nucleation during changes in the angle cp from 0° to 30°. The Thompson Triangle 17 shown in Fig. 
4 is now used to illustrate the transition from twinning to full dislocation nucleation in the active 
(ll l) slip plane. At cp = 0° (Fig. 4a) the direction of the tensile load, the y-direction, is parallel to 

the twinning direction [112]. In this case, both nucleated partial dislocations, c/1 and c/2, are of 
the same type with Burgers vectors in the [112] direction. This is the situation when a twin is 
formed. At cp = 13.90° (Fig. 4b), the Thompson triangle is slightly rotated in a clockwise 
direction, so that the y-direction of the tensile load is slightly inclined to the twin direction [112], 
but the inclination angle cp is too small to affect the dislocation nucleation process. The 
dislocation nucleation process is similar to that in Fig. 4(a), and twinning still takes place. 
However, further increase of the angle cp does change the deformation mode from twinning to 
full dislocation emission. At cp = 16.10° (Fig. 4c), the second dislocation d2 is nucleated in the 
[l2 1 ] direction. Both partial dislocations, d\ and c/2, combined to form a full dislocation c/3 with 

a resultant Burgers vector along the [011] direction. The same situation occurs for the remainder 
of the cp angles as indicated in Figs. 4(d-f). 

A simplified but illustrative explanation of the change in the deformation mode from twinning 
to full dislocation emission can be given. Deformation twinning produces displacements in the 
crystal lattice along the <1 12> family of directions through the consecutive emission of parallel 
Shockley partial dislocations 15 such as c/1 and c/2 in Fig. 4(a) and Fig. 4(b). Full dislocation 
emission produces displacement along the <11 0> family of directions by the emission of two 
Shockley partial dislocations in non-parallel <112> directions with a combined Burgers vector 
along one of the <11 0> directions. Schematically, this is shown in Figs. 4(c-f) where the 
combined dislocation c/3 is formed as a sum of the two non-parallel partial dislocations, c/1 and 
c/2. When cp < 15°, the tensile y-direction is closer to the [112] direction (Fig. 4a and Fig. 4b) and 
the displacement takes place in the [112] direction, which results in twinning. When cp > 15°, the 
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tensile v-direction is closer to the [Oil] direction (Figs. 4c-f), and it is more favorable for the 
displacement to take place in the [Oil] direction. In order for this to happen, the second partial 
dislocation d2 is nucleated in the [l2 l] direction so that the total Burgers vector of d\ and d2 
produces a displacement along the [Oil] direction. 



Figure 4: Crystallographic orientations of the (ill) slip system at the crack tip in Thompson schematics for the set 
of orientations given in Table 1. 


A more rigorous analysis of the process is given by Tadmor and Hai 3 . They calculate the 
critical stress intensities, K^ m and Kf ( d , required to nucleate a pair of parallel (for twinning) or 
non-parallel (for full dislocation emission) Shockley partials under a given direction of loading 
and crystal orientation. If K^ m < Kf ( f then twinning is the more energetically favorable process 
and is more likely to occur. If K^ 1 < then full dislocation emission is the energetically 


11 

American Institute of Aeronautics and Astronautics 


more favorable result. Based on this consideration, Tadmor and Hai 3 proposed a criterion for 
twinning against full dislocation emission as the ratio between K^ 1 and K 


For the case of pure mode I loading, the calculations show 3 that, while K']’f and K™" 
separately depend on the tilt angle 0 (Fig. 2), their ratio T(cp ) does not. In this way, the results for 
0 = 90°, presented in this work, are representative for all tilt angles scenarios. Using the 
methodology by Tadmor and Hai 3 to calculate Kf^{cp) and K'" c m (cp ) for the energetically most 
favorable pair of partial dislocations, the dependence of their ratio T(cp ) on cp for aluminum is 
given in Fig. 5. The material parameters that enter into the calculation 3 are the stable and 
unstable stacking fault energies, written as y s f and y us , respectively, as defined by Rice 18 , and the 
unstable twinning energy, y ut , as defined by Tadmor and Hai 3 . These energies for the curve, 
presented in Fig. 5, were derived from the interatomic potential 9 at T = 0 K. The specific values 
are given as 9 : y s f = 0-146 J/m 2 , y U s = 0-168 J/m 2 , and 4 y ut = 0.21 J/m 2 . It is expected that their 
values at T = 300 K are slightly lower, and the curve on Fig. 5 may not be quantitatively 
representative for T = 300 K, but it is given here for qualitative analysis. 

From a qualitative point of view, Fig. 5 shows a maximum value of T(cp ) at cp = 0° and a 
minimum at cp = 30°. Thus, twinning is expected at cp = 0° and full dislocation emission is 
expected at cp = 30°, which is in agreement with the simulation results in Fig. 4(a) and Fig. 4(f). It 
should be noted however, that the range where T(cp)> 1 in Fig. 5 is very narrow (about ±1°), 
while the MD simulations show a predominant twinning in a broader range up to cp = 13.9° (Fig. 
4b). One possible explanation for this quantitative disagreement could be that the analytical 
curve in Fig. 5 is given for T = 0 K, while the simulations were performed for T = 300 K. A 
more likely explanation is that the analytical curve is based on a minimum energy consideration 
and represents a static case, while the simulations are dynamic and were performed at a very 
high strain rate of 5xl0 7 s' 1 . The mechanisms by which the strain rate may affect the twin- 
dislocation transition, as found in the atomistic simulations, are discussed in the next subsection. 



Rdisi j > ] twinning 

Kj " c in T < 1 dislocation 


12 

American Institute of Aeronautics and Astronautics 



1.1 



III.2. Twin-dislocation transition with stress intensity 

The effect of the high strain rate, inherently present in molecular-dynamics simulations, on 
the crack tip twinning in aluminum has recently been discussed by Warner et al. 5 The main 
consideration in that work 5 , based on transition state theory 19 , was that to observe a dislocation 
nucleation from the crack tip in an MD simulation time frame (a few nanoseconds), the 
dislocation nucleation rate must be on the order of 10 9 -10 10 events per second. To achieve this 
high nucleation rate, much larger stress intensities have to be applied as compared to 
experimental conditions, where a nucleation rate of 1 0°- 1 0 1 is usually enough to observe 
dislocation nucleations in the time frame of seconds. The role of the stress intensity on the 
nucleation rate is accounted for in the Rice and Beltz model 20 of dislocation nucleation at a crack 
tip. The analysis 5 has shown that a transition stress intensity for twinning Kn exists, at which the 
deformation mode changes from full dislocation emission to deformation twinning. For the case 
cp = 0°, considered as a typical twinning case 5 , K n was found to be substantially lower than 
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values of K\ that are commonly observed in MD simulations of crack propagation. For a full 
dislocation emission to take place, the condition K\ < K\j has to be fulfilled. However, according 
to transition state theory 5 ' 19 , the dislocation nucleation rate at K\ < K\j is so low that to record a 
dislocation emission event the observation time in the MD simulation has to be increased by 
more than four orders of magnitude to the level of 1CT 5 s (equivalent to a strain rate of 10" 3 s’ 1 for 
1% deformation). To achieve this slowdown, an advanced coupled multiscale-accelerated 
dynamics technique, which involved scaling on both, space and time scales, was applied 5 . 

Warner et al. calculations 5 for the transition from deformation twinning to full dislocation 
nucleation were found to agree qualitatively with the theoretical predictions, but substantial 
quantitative differences were also apparent, indicating that further study is needed to fully 
understand the process. The complexity of the applied multiscale-accelerated dynamics 
technique by Warner et al. 5 made its use more difficult in a more systematic study of the 
problem. A different approach is presented in the present paper, which allows the twin- 
dislocation transition to be reproduced and studied at the MD-accessible strain rate of 5x1 0 7 s’ 1 . 
This permits the more computationally affordable conventional MD technique to be used and 
avoids some possible artifacts of applying an accelerated-dynamics technique. 

The new approach is based on the observed MD twin-dislocation transition with the change 
of crystal orientation as described in Sec. III. 1 . The observed transition from twinning to full 
dislocation emission shown in Fig. 4 can be interpreted as a result of increasing the transition 
stress intensity, Ku, as cp varies from 0 to 30°. Thus, at some intermediate angle, cp, between 
13.9° and 16.1°, K\j becomes higher than the applied stress intensity A), and twinning is replaced 
by full dislocation emission. If this is indeed the case, then increasing the applied stress intensity 
above Kn would reinitiate the twinning mechanism, whereas full dislocation emissions were 
initially observed at lower stress intensities. This scenario is easier to test by MD simulations 
than the case studied by Warner et al. 5 , because in an MD simulation, an increase of K\ is much 
easier to achieve than a decrease, as it avoids the use of accelerated dynamics techniques. The 
drawback is that increasing the stress intensity may induce additional deformation modes, such 
as the activation of a secondary slip system, as seen in Fig. 6, which will be discussed next. 

To examine the dependence of the twin-dislocation transition on increasing stress intensity, 
the first two values of cp, equal to cp = 16.10° and cp = 19.11°, were selected because they emitted 
full dislocations. The systems were subject to four levels of uniaxial tensile strain of e yy = 0.7%, 
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0.8%, 0.9%, and 1.0% corresponding to K\ = 0.34, 0.39, 0.44, and 0.49 MPa-m 1/2 , respectively. 

1 /9 

The results for cp = 16.10° are shown in Fig. 6. At K\ = 0.34 MPa-m (Fig. 6a), the crack emitted 

i j'y 

a full dislocation and at K\ = 0.39 MPa-m , a twin was nucleated (Fig. 6b). Twin nucleation 
continued to be systematically observed at higher loads of K\ = 0.44, and 0.49 MPa-m (Figs. 6c 
and 6d). From these results, it follows that the transition stress intensity at cp = 16.10° is 
approximately Ky r = 0.39 MPa-m ' . Similar behavior is found for cp = 19.11°, as shown in Fig. 7, 
but twinning starts at a higher load of K\ = Kn = 0.49 MPa-m . Thus, an increase in cp results in 
an increase in Kn by the same mechanisms that were observed in Fig. 6. The high stress at which 
twinning starts, induces the emission of partial dislocations on a secondary slip plane as seen in 
Figs. 6(c,d) and Figs. 7(c,d). 



Figure 6: Simulation snapshots of the initiation of plastic processes at the crack tip for cp = 16.10° at different stress 
intensities K u as indicated in the figure. 
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The presented results confirm the conclusion of Warner et al. 5 that a transition stress intensity 
Kyi exists, below which the crack emits full dislocations, and above which deformation twinning 
takes place. The results also show that Kn is a function of the twist angle (p. At (p = 0, Kn is too 
low for a dislocation nucleation event to take place at typical MD time scales and full dislocation 
emission cannot be seen in a conventional MD simulation at this orientation. At orientations for 
which Kn is above the simulated K\, full dislocation emission becomes the dominant crack 
propagation mode. When K\ is increased above Kn, deformation twinning is again observed. 



Figure 7: Simulation snapshots of the initiation of plastic processes at the crack tip for cp = 19.1 1° at different stress 
intensities K u as indicated in the figure. 


IV. Summary 
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The simulation method for this study is based on the embedded statistical coupling method 
(ESCM). The ESCM is a MD-FEM coupling method, which uses statistical quantities such as 
statistically averaged atomic displacements to connect the atomistic MD domain with the 
continuum FE domain. As a consequence, the ESCM technique makes possible relatively large 
atomistic domains (10 5 to 10 7 atoms) to be embedded in a continuum domain of micrometer 
dimensions. This capability is essential for the present study because it minimizes the finite size 
effects by allowing edge cracks of 0.45 pm size to be modeled while preserving the atomic 
representation of the material at the crack tip. The ESCM method so far has been successfully 
applied in the case of intergranular fracture, where the crack propagation direction was 
predefined by the existing grain boundary. The presented study is the first application of the 
ESCM approach for the case of transgranular fracture. 

Dislocation processes at a crack tip in a single crystal of pure aluminum were studied. Two 
types of processes were found to take place: deformation twinning and emission of full 
dislocations (i.e. dislocation slip at the crack tip). Studying the transition between these two 
processes reveals the existence of a transition stress intensity, Kn, below which the crack emits 
full dislocations and above which deformation twinning becomes dominant. The transition stress 
intensity depends on the crystallographic orientation of the crack front in mode I loading, 
through the twist angle cp between the tensile direction and the twinning [112] crystallographic 
axis. A minimum value of Kn is reached at cp = 0° where twinning becomes the dominant crack 
propagation mode in MD simulations due to their inherently very high strain rates of 10 7 s' 1 or 
higher. A maximum value of Kn at cp = 30° defines the region of full dislocation emission at 
typical MD time scales. To be consistent with experimental observations, where deformation 
twinning at the crack tip in aluminum is rarely observed, this study suggests that crystallographic 
orientations close to cp = 30° should be used for atomistic characterization of crack tip plastic 
processes in aluminum. Orientations close to cp = 0° should be avoided as they produce the 
artifact of deformation twinning, that can substantially alter the predicted fracture parameters, 
such as the peak stress of debonding and energy of decohesion. 
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